use strict;
use warnings;
use Data::Dumper;
use Test::More 'no_plan';
use File::Temp qw(tempdir);
use Cwd;
use Bio::Rambam::Proteomics::old::Quant;
#use Bio::Rambam::Proteomics::old::Formats::kkMzML;
#use Bio::Rambam::Proteomics::old::Formats::kkMzMLSpectrum;
#use Bio::Rambam::Proteomics::old::Formats::OmssaFile;
#use Bio::Rambam::Proteomics::Experiment;

#my $original = cwd;
my $psms_iniciales=0;
my $borrados_score=0;
my $borrados_son_MS1=0;
my $borrados_no_quant=0;
my $siguen_vivos=0;

my $search_file="/home/david/projects/jaen/Santos_TMTGrupo1_22102013.mgf.txt";
my $mzMLfile="/home/david/projects/jaen/Santos_TMTGrupo1_22102013.mzML";
my @csv_files;
my $omssa_data_file="./testdata/Santos_TMTGrupo1_22102013.mgf.txt";

my $t1=time();
my ($i126,$i127,$i128,$i129,$i130,$i131);
($i126,$i127,$i128,$i129,$i130,$i131)=(0,0,0,0,0,0);

push @csv_files,$omssa_data_file;
my $experiment=Bio::Rambam::Proteomics::Experiment->new(name=>'Data 1 and 2 corrected');
$experiment->process_omssa_csv_files(\@csv_files);
$experiment->process_quant_info(method=>'TMT6plex',mzml_file=>$mzMLfile,ident_scan=>'CID',quant_scan=>'HCD',quant_to_ident=>+3);
$experiment->assign_unfiltered_psms_to_proteins();

my %proteins=$experiment->get_proteins;
my %psms=$experiment->get_psms;
foreach my $prot(keys %proteins){
	my $used_psms=0;
	my (@r1,@r2,@r3,@r4,@r5);
	my $protein=$proteins{$prot};
	next unless $protein->num_psms_ref>2;
	next unless defined $protein->quant_info;
	#print $prot,"\t->\t".$protein->num_psms_ref."\t\n";
	foreach my $qik(keys %{$protein->quant_info}){
		my $psmquantinfo=${$protein->quant_info}{$qik};
		my %quants=%{$psmquantinfo->quantInfo};
		$i126=$i126+$quants{126};
		$i127=$i127+$quants{127};
		$i128=$i128+$quants{128};
		$i129=$i129+$quants{129};
		$i130=$i130+$quants{130};
		$i131=$i131+$quants{131};
		if($quants{126}==0){
			#print "----\n";
			next;
		}
		$used_psms++;
		my $ratio_1=sprintf("%.3f",($quants{127}/$quants{126}));
		my $ratio_2=sprintf("%.3f",($quants{128}/$quants{126}));
		my $ratio_3=sprintf("%.3f",($quants{129}/$quants{126}));
		my $ratio_4=sprintf("%.3f",($quants{130}/$quants{126}));
		my $ratio_5=sprintf("%.3f",($quants{131}/$quants{126}));
		my $psm=$psms{$qik};
		print $qik."\t";
		print $prot."\t".$psm->peptide->sequence."\t";
		print $ratio_1."\t".$ratio_2."\t".$ratio_3."\t".$ratio_4."\t".$ratio_5."\n";
		push @r1,$ratio_1;
		push @r2,$ratio_2;
		push @r3,$ratio_3;
		push @r4,$ratio_4;
		push @r5,$ratio_5;
		#print Dumper $psmquantinfo;exit;
		#print ${{$protein->quant_info}{$qik}->quantInfo{128};
	}
	next unless $used_psms>2;
	my $total_ratio_1=sprintf("%.3f",average(\@r1));
	my $total_ratio_2=sprintf("%.3f",average(\@r2));
	my $total_ratio_3=sprintf("%.3f",average(\@r3));
	my $total_ratio_4=sprintf("%.3f",average(\@r4));
	my $total_ratio_5=sprintf("%.3f",average(\@r5));
	my $cv1="NA";
	my $cv2="NA";
	my $cv3="NA";
	my $cv4="NA";
	my $cv5="NA";
	$cv1=sprintf("%.3f",100*(stdev(\@r1)/$total_ratio_1)) unless $total_ratio_1==0;
	$cv2=sprintf("%.3f",100*(stdev(\@r2)/$total_ratio_2)) unless $total_ratio_2==0;
	$cv3=sprintf("%.3f",100*(stdev(\@r3)/$total_ratio_3)) unless $total_ratio_3==0;
	$cv4=sprintf("%.3f",100*(stdev(\@r4)/$total_ratio_4)) unless $total_ratio_4==0;
	$cv5=sprintf("%.3f",100*(stdev(\@r5)/$total_ratio_5)) unless $total_ratio_5==0;
#	print $prot."\t".$protein->name.$protein->num_psms_ref."\t$used_psms\t";
#	print $total_ratio_1."\t".$cv1."\t";
#	print $total_ratio_2."\t".$cv2."\t";
#	print $total_ratio_3."\t".$cv3."\t";
#	print $total_ratio_4."\t".$cv4."\t";
#	print $total_ratio_5."\t".$cv5;
#	print "\n"
}


#print "$i126\n";
#print "$i127\n";
#print "$i128\n";
#print "$i129\n";
#print "$i130\n";
#print "$i131\n";



#print Dumper $experiment;


exit;
#$ave = &average(\@array);
##$ave = Math::NumberCruncher::Average(\@array); //OLD WAY
#$std = &stdev(\@array);
##$std = Math::NumberCruncher::StandardDeviation(\@array); //OLD WAY
sub average{
        my($data) = @_;
        if (not @$data) {
                die("Empty array\n");
        }
        my $total = 0;
        foreach (@$data) {
                $total += $_;
        }
        my $average = $total / @$data;
        return $average;
}
sub stdev{
        my($data) = @_;
        if(@$data == 1){
                return 0;
        }
        my $average = &average($data);
        my $sqtotal = 0;
        foreach(@$data) {
                $sqtotal += ($average-$_) ** 2;
        }
        my $std = ($sqtotal / (@$data-1)) ** 0.5;
        return $std;
}


#my %proteins=$experiment->get_unfiltered_proteins();




########    Ahora vendria esto de abajo:
#$experiment->process_itraq_info();
#my @peptides=$experiment->get_peptides();
#my %proteins=$experiment->get_proteins();

#tienes que tocar, trayendolo de la libreria Proteomics:
#experiment->process_itraq_info
#$protein->set_itraq_stats($R);
#(...)
#y tira millas, que no te pase nada....




exit;
















#END {
#﻿  chdir $original;
#}